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ABSTRACT 

We present a new method to find the transfer function (TF) of the broad-line region 
in active galactic nuclei. The subtractive optimally localized averages (SOLA) method 
is a modified version of the Backus-Gilbert method and is presented as an alterna- 
tive to the more often used maximum-entropy method for the inversion of variability 
data. The SOLA method has been developed for use in helioseismology. It has been 
applied to the solar oscillation frequency splitting data currently available to deduce 
the internal rotation rate of the sun. The original SOLA method is reformulated in 
the present paper to cope with the slightly different problem of inverting time series. 
We use simulations to test the viability of the method and apply the SOLA method 
to the real data of the Seyfert-1 galaxy NGC 5548. We find similar TFs for these data 
as previous studies using the maximum-entropy method. We confirm thereby previ- 
ous results while simultaneously presenting an alternative and independent inversion 
method. Moreover, we do not find significant negative responses in the TF. The inte- 
gral of the TF, an important quantity measuring the total observed line-processing in 
the broad-line region, is correctly reproduced by the SOLA method with high accu- 
racy. We investigate the effects of measurement errors and how the resolution of the 
TF critically depends upon both the sampling rate and the photometric accuracy of 
the data. 
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REVERBERATION MAPPING 



Mathematically, the concept of reverberation mapping 
reduces to the inversion problem 



The broad-Une region (BLR) in active galactic nuclei 
(AGNs) is too small to be resolved spatially with even the 
largest planned telescopes. Indirect methods must therefore 
be employed to obtain information about its structure and 
dynamics. One way of doing so is via reverberation or echo 
mapping, a technique possible for variable sources only. 



where L{v,t) is the observed emission-line light curve, v the 
projected velocity with v = the line centre, C{t) the ob- 
served continuum light curve, and 4'(u,r) the transfer func- 
tion (TF) of the BLR. The TF is the quantity that holds 
the information about the geometry, kinematics and physics 
of the BLR and is thus of central importance. To obtain 
higher S/N in the time series, Eq. (1) is usually integrated 
over projected velocity v which reduces it to the total-flux 




(1) 



The broad emission lines originating in the BLR are 
photoionized by a small central continuum source (see Net- 
zer 1991 for a detailed review of the standard model). In 
many AGN, especially the low-luminosity ones, this contin- 
uum source is observed to be variable in its flux. The broad 
emission lines respond to the continuum variations, albeit 
with a time delay of several days due to light-travel-time 
effects through the BLR. The high velocities of the clouds 
comprising the BLR redistribute line photons in wavelength, 
resulting in the broad emission-line profiles. Therefore, the 
combination of flux and proflle variations of the emission 
lines in response to the ionizing-continuum variations can 
be used to map the phase space of the BLR (Blandford & 
McKee 1982; see Peterson 1993 and references therein for 
an extensive review). 



case 




To date two inversion methods have been applied to 
AGN variability data to invert Eq. (2) and so deduce the 
TF. These are the Fourier-transform method (FTM) (e.g., 
Bracewell 1986; Maoz et al. 1991) and the maximum-entropy 
method (MEM) (SkiUing & Bryan 1984; Home, Welsh & 
Peterson 1991; Krolik et al. 1991). The FTM is very noise- 



(2) 
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sensitive and cannot handle irregularly sampled light curves 
properly. Only for the galaxy NGC4151 has a TF been de- 
rived via this method (Maoz et al. 1991) but without a clear 
result. The MEM on the other hand has boon successful in 
finding a number of TFs for several AGN: NGC 5548 (Horne, 
Welsh & Peterson 1991; Krolik ct al. 1991; Peterson 1993; 
Peterson et al. 1994), Mkn590 (Peterson et al. 1993) and 
NGC 3516 (Wanders & Horne 1994). A third method, based 
on a linear rcgularization method (cf. Press et al. 1992) is 
currently under development (Krolik 1994). 

In this paper we will outline another method for the in- 
version of Eq. (2), not hitherto applied to AGN. This method 
is a modified version of the Backus-Gilbert method (Backus 
& Gilbert 1967, 1968, 1970) in use in gcoscismology. More 
recently the method has been applied to helioseismology, 
specifically to the solar-oscillation frequency-splitting data 
to deduce the internal rotation rate of the sun. A comparison 
of various methods in use in this field, including the Backus- 
Gilbert method, can bo found in papers by Gough (1985) 
and by Christensen-Dalsgaard, Schou & Thompson (1990). 
The Backus- Gilbert method is also sometimes referred to as 
optimally localized averages (OLA) or multiplicative OLA 
(MOLA). The MOLA methods are generally considered not 
to be cost-effective since they require numerous matrix in- 
versions of large matrices. This led Pijpers & Thompson 
(1992, 1994) to reformulate the OLA method into a fast 
method called the subtractive optimally localized averages 
(SOLA) method. In Sect. 2 this method will be outlined in 
more detail and the transformation from the helioseismolog- 
ical formulation to the reverberation-mapping concept will 
be made. Section 3 describes on which points the application 
to reverberation mapping differs from helioseismology, and 
problems that are specific to AGN time-series analysis will 
be addressed. Section 4 discusses simulations done on real- 
istic (i.e. noisy and irregularly sampled) artificial data and 
shows that recovery of the TF by the SOLA method can be 
done with currently available data. In Sect. 5 wo then apply 
the SOLA method to the real data of the Seyfert-1 galajcy 
NGC 5548 and present its TF for four years of monitoring 
data. We show that the presently derived TFs are consis- 
tent with the previously derived ones using the MEM and 
the observed global changes in the TF from year to year arc 
reproduced by the SOLA method. These global changes are 
both changes in the position of the peak of the TF and in the 
integral of the TF. The evolution of the BLR in NGC 5548 
on the dynamical time scale is thereby confirmed. In Sect. 6 
we summarize our results. 



2 THE SOLA METHOD 

The strategy of the SOLA method is to find a set of linear 
coefficients which, when combined with the data, produce 
the unknown convolved function under the integral sign. In 
the application to reverberation mapping these coefficients 
will be multiplied with the individual line fluxes in the mea- 
sured time series and then summed to obtain an estimate of 
the TF at a certain lag tq ■ 

Consider again Eq. (2) but now with the line fluxes 
explicitly written as a time series consisting of A'' individual 
measurements L{ti) with an associated measurement error 



SL{ti) 



L{ti) - 6L{t 



dr *(T)C(ti - t) \fi=l,...,N. (3) 



The lower limit of the integration is set to which implies 
that the line flux is assumed to respond to the continuum 
variations and not vice versa. If the contiimum flux C would 
be known for all values of t and would be error-free this 
problem would be equivalent to the problem of inverting 
helioseismology data. In that case an averaging kernel K, 
would be constructed : 



^(T-0,r) = Y^Ci{To)C{ti-T) 



(4) 



which is strongly peaked around ro and which is normal- 
ized : 



oo 

/ 



dr K,{to,t) = 1. 



(5) 



An estimate *(to) of the TF at time lag r = To is then 
obtained by combining the line fluxes with this same set of 
coefficients 



$(ro) = Y,Ci{To)L{U) = Ut K:(ro,T)*(r) 



(6) 



+ ^Ci{To)SL{U). 



The first term on the right of Eq. (6) is a weighted average 
of ^ in which IC is the weighting function. The other term 
expresses the propagation of the errors in the line fluxes. In 
the SOLA method the aim is to make the kernel K, resemble 
some chosen target form T while at the same time moderat- 
ing the effect of the errors in the line fluxes. This is achieved 
by minimizing 



oo 

j At [1C{to,t) -T{TQ,T)f + n'y^ EijCiCj 



(7) 



subject to the constraint (5). Here E is the error variance- 
covariancc matrix of the observed line fluxes, /i is a trade- 
off' parameter which may be chosen according to the relative 
desirability of maJiing the first and second terms in (7) small. 
The ideal target function for the averaging kernel would be 
a delta function but it is clear that an infinite resolution 
cannot be attained in practice. A convenient form of a target 
function with an adjustable width is 



T = 



1 
/A 



exp 



(8) 



Here / is a normalization factor which is introduced to make 
the total integral of T equal to unity. If A -C 1 the factor / 
reduces to a/tF- 

The minimization leads to the matrix equation (Pijpers 
& Thompson 1994) 



A • c(ro) 



.(ro). 



(9) 
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Here c(to) is the vector of linear coefficients which is un- 
known, with an extra (A^^-l- l)th element which is a Lagrange 
multiplier. The elements of the symmetric matrix A are : 

dT C{U - T)C{tj - t) + tiEij {i, j < N) 

l^drCiU-r) {i<N,j = N+l) 
dr - r) (i<iV,i = iV + l) 
{i=j = N+l) 

The elements of the vector v are given by 



(11) 



, , . / dr C{ti-T)T{ro,r) {i < N) 

Vi{To) = J 

1 {i = N + l) 

Note that the elements of the matrix A do not depend on 

the point tq at which an estimate of the TF will be deter- 
mined. The matrix has to be inverted only once. After this 
the different vectors v(ro) can be combined with the inverse 
matrix with relatively little computational effort, to scan the 
TF over the range of To of interest. It should also be noted 
that if a different error weighting /i is chosen the matrix 
elements of A can be calculated quickly but a new matrix 
inversion is required. 

In practice one finds that there is a trade-off between 
the error weighting /x and the resolution width A. Increasing 
/i produces a decreasing error in the TF propagated from 
the line flux measurements but the averaging kernel IC will 
depart more and more from the target form T unless A is 
also increased. This departure of IC from the prescribed form 
introduces a 'systematic error' in the determination of the 
TF which is undesirable. A more detailed discussion of this 
trade-off can be found in the paper of Fijpers & Thompson 
(1994). 



3 ASPECTS OF TIME-SERIES ANALYSIS 
WITH THE SOLA METHOD 

3.1 Adapting SOLA 

It is clear that for the problem of time series analysis dis- 
cussed here the original formulation of SOLA cannot work 
since the continuum fluxes are known only for discrete points 
in a finite time interval. An extra problem is that this time 
series for the continuum flux also has measurement errors. 
These problems will be dealt with in several steps leading 
to a working algorithm for SOLA inversion of time series. 

3.2 The finite extent of the continuum time series 

In Fig. 1 an example is shown of a time series of mea- 
surements of the continuum flux. The index i counts the 
measurements. For each added measurement the whole time 
series is replotted with an arbitrary vertical offset. The time 
series is plotted as a function of delay time t so that the 
most recent measurement is always the left-most point at 
r = 0. The dashed lines represent the unknown part of the 
light curve before the first measurement. For this time series 
of 10 measurements the 10 curves shown in Fig. 1 would be 
the set of base functions C that appear in Eqs. (10) and (11). 
Ftom Fig. 1 it is clear that it is not possible to evaluate the 
integrals since there is always a part of the base functions 
C which is not known since it falls before the first mea- 
surement. The solution to this problem is to decide before 



figure 1 Example of a time series with 10 measurements. The 
index i counts the consecutive measurements. The whole time 
series is re-plotted as a function of delay time r for each added 
measurement, with an arbitrary vertical offset. The dashed lines 
represent the flux before the first measurement in the series 

starting the inversion to put the upper limit of the integra- 
tions to a certain finite time lag Tmax- This implies that the 
TF will be assumed to be negligible for time lags t > Tmax- 
In essence this is an assumption regarding the maximum size 
of the BLR. If this parameter is unknown one must exper- 
iment with various values for Tmax and re-do the inversion 
for each value. 

The maximum time lag Tmax is indicated in Fig. 1 by 
the vertical dash-dotted line. Even with a finite value of 
Tmax there will still be some base functions that are known 
for only part of the range [0, Tmax] so these will have to be 
dropped completely. This reduces the number of indepen- 
dent base functions. It is important not to choose Tmax too 
large because then the number of independent base func- 
tions becomes very small which has an adverse effect on the 
error in the reconstructed TF and also on the resolution that 
can be attained. 

It is also important not to choose Tmax too small since 
then each base function is determined by only a very few 
points. This means that the integrals in (10) and (11), which 
are calculated on the basis of an interpolation, become un- 
certain. Note that the discarding of some of the base func- 
tions does not mean that the first continuum measurements 
in the series are discarded. These are all used but shifted 
to the right in the base functions for higher i. However, if 
it is assumed that the measurements of continuum flux and 
line flux are always simultaneous there is some data that is 
obtained observationally which cannot be used and must be 
discarded. If the first M base functions are discarded the 
first M line-flux measurements are also discarded. Since the 
line flux of these first M measurements (paxtly) responds 
to variations in a part of the continuum flux that has not 
been measured, these measurements are of limited use for 
the inversion. Because the TF is only calculated at one tq 
at a time, one could, in principle, use a floating Tmax, i.e., 
Tmax is short if To is short, such that fewer emission- line data 
are ignored. However, this would mean the matrix inversion 
(9) must be done for each different Tmax which is expensive. 
Furthermore, with finite resolution, boundary effects play a 
significant role when Tmax becomes of the order of twice the 
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resolution or less. All in all, it is not advantageous to do this 
and in this paper we will keep Tmax fixed and drop the first 
M line-flux measurements. 

3.3 The discrete sampling of the continuum time 
series 

Even after limiting the inversion to a finite interval in r 
the base fiinctions arc not completely determined. The con- 
tinuum flux is only known on a discrete set of points. The 
second step towards a working algorithm is to interpolate 
between the individual measurements. It is important to re- 
alize that this has an impact on the resolution that can be 
attained in the determination of the TF. It implies that there 
is a minimum to the width A of the target function (8) even 
with error-free data. 

At this point it is useful to discuss the 'intrinsic' res- 
olution that could be attained with a series of error-free 
continuum-flux measurements. If the sampling rate of the 
continuum flux is strictly regular the base functions also are 
sampled strictly regularly. It is immediately clear that the 
minimum possible resolution is then of the order of the time 
interval between samplings. If the sampling is irregular the 
result is less obvious. It is useful to distinguish at this point 
between a regular time series with 'missing data' and a fully 
irregular or non-redundant time series. 

In a time scries of the first type the sampling can be 
thought of as taken on integer multiples of some minimum 
time interval Aimin- The 'missing data' are thou all those 
integer multiples of At for which a measurement is not avail- 
able. In this case in the diagram of Fig. 1 a set of K regularly 
spaced vertical lines can be drawn with K = Tmax/Atmin- 
All measured points of all the base functions will fall on 
those K lines. Considered in this way the minimum resolu- 
tion width A with this type of irregular sampling must be 
of the order of Atmin. If a fraction / of the data is missing 
in a time scries with this imaginary sampling rate Atmin the 
length of the total time series must be larger by a factor 1/f 
than a similar time series without missing data, in order to 
really attain something close to this resolution. 

In a non-redundant time series there is no minimum 
sampling interval for which all measurements fall on integer 
multiples of that time interval. Note that this does not re- 
quire that some samples are taken with an infinitely small 
time interval between them. Consider again Fig. 1, where 
now the interval between and Tmax is rescaled so that it is 
mapped onto the interval [0, 1]. In the limit of an infinitely 
long non-redundant time series the set of r/rmax for which 
at least one measurement is available would lie dense in the 
unit interval. In other words the distribution of these points 
would have the same properties as the rational numbers on 
the unit interval. For an error free time series of this type 
the minimum resolution width is inversely proportional to 
the total length of the time series. 

The time scries of measurements of fluxes of variable 
AGN as they are usually reported fall into the flrst category 
of irregularly sampled data. The minimum time interval here 
is 1 day. However, one should not expect to attain this lim- 
iting resolution of 1 day for two reasons. The first reason 
is the finite error in the measurements of the line fiux. This 
leads to the usual trade off between errors in the TF and res- 
olution of the TF, mentioned in the previous section. The 



second reason is finite error in the continuum fiux measure- 
ments. These errors give rise to an error or uncertainty in 
the integration kernel JC that is constructed. This error in- 
creases for a decreasing width of the target kernel T which 
effectively limits the resolution width to a minimum which 
may lie well above the minimum of 1 day for error-free data. 

There is one other reason why this limiting resolution 
may not be attained which has nothing to do with the sam- 
pling strategy or the measurement errors. It is clear that if 
the continuum does not vary, then it is impossible to deduce 
anything about a TF. As an extension to this it is not hard 
to see that if there is a finite minimum time scale on which 
the continuum source varies, which is larger than the sam- 
pling time interval, one should not expect to attain a better 
time resolution in the TF than that typical time scale of the 
continuum source. This should in principle be reflected in 
an impossibility for the algorithm to find a good kernel K, 
with a smaller width than this time scale (cf . Fig. 7) . 

3.4 Interpolation of the continuum time series 

The practical problem left is to find an interpolation scheme 
for the continuum fiux time series so that the integrals of 
(10) and (11) can be calculated. This interpolation must 
not introduce spurious variations due to the errors in the 
continuum flux. When tested on simulated time series a cu- 
bic spline was found to introduce spurious oscillations with 
amplitudes of up to four times the formal errors on the indi- 
vidual measurements. Such spurious structure in the contin- 
uum flux would of course have no counterpart in the line-flux 
measurements. In this way an artiflcial de-correlation would 
be found for some determinations of the TF, effectively forc- 
ing the TF towards values close to 0. This would happen in 
particular for values of to for which none or only a few base 
functions have a measured point. 

The interpolation scheme decided upon here is a 
Savitzky-Golay smoothing filter for the contirmum time se- 
ries, adapted to an irregular sampling interval (cf. Press et 
al., 1992). With this algorithm a least-squares fit of a low or- 
der polynomial to a moving window of 2A's -I- 1 points in the 
time series is constructed. Except at the beginning and end 
of the entire series these points are always distributed sym- 
metrically around the point for which the polynomial will 
be used in the interpolation. An interpolating polynomial 
constructed for the continuum light curve around time ti is 
then used only in the interval [{ti-i+ti)/2, {ti + ti+i)/2]. For 
the A's points at both ends of the time series the number of 
points used is reduced one at a time until only A's -|- 1 points 
are used to construct an interpolating polynomial for the 
first and last point in the time series. In this case the order 
of the polynomial is also reduced to make sure that the error 
in the coefficients of the polynomial remains under control. 
A choice of a moving window of 5 points (i.e. Ns = 2) and a 
quadratic polynomial are found to perform satisfactorily. At 
the boundaries the polynomial was reduced to linear. The 
choice of the number of points in the smoothing window and 
the degree of the polynomial iu general must depend on the 
ratio of the real point-to-point variations in the continuum 
light curve and the statistical measurement noise. 

Note that this interpolation scheme is not continuous 
for all T in the range [0, Tmax], contrary to a cubic spline. 
This is of no importance for the algorithm but it can show 
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up in the reconstructed kernel IC which then appears some- 
what jagged. Note also that the polynomial is intended to 
produce an interpolation without spurious oscillations. It 
is not intended to significantly reduce the errors from the 
data, for which a much larger window should be chosen. 
Suppressing the errors from the continuum is of secondaxy 
importance since this can better be done by increasing the 
total length of the time series, and therefore the number of 
base functions, which only requires additional observations. 
The cost of smoothing with a window with more points is 
a degraded limiting resolution unless the sampling is done 
at a higher rate. It is certainly impossible to get a higher 
sampling rate in the part of a time series that is already 
obtained. Also there arc practical limitations to the highest 
sampling rate in observing campaigns, whereas extending a 
time series is always possible in principle. 

The resolution limit imposed by this scheme is not uni- 
form since it is proportional to the length in time of the 
smoothing window. A guide-line for the smallest value of 
A is found by equating the FWHM of the Gaussian target 
function with the width of the window: 

Amin ~ ^^°jL^ Atav > Atniin (12) 

2Vln2 

Here Atav is the average time interval between individual 
flux measurements. As an example consider a time scries of 
100 measurements spread over one year. The limiting value 
for the resolution width would be A ~ 11 d. 

3.5 Error propagation 

Combining the steps in the previous sections leads to a work- 
ing inversion algorithm. Such an algorithm is not much use 
without an estimate of the error in the determination of 
the TF, which are simply propagated measurement errors. 
In the original SOLA formulation the error propagation is 
relatively simple. The TF is obtained from linear combina- 
tions of the line fluxes and therefore the variance of the TF 
is simply 

a'(*(To)) = c^(to)-El-c(to) (13) 

Here El is the error variance-covariance matrix of the line- 
flux measurements. Usually the errors in the individual line- 
flux mcEisurements are assumed to be uncorrelated which 
means that the matrix El is diagonal. 

In the new formulation of SOLA the coefficients c 
should also have an error associated with them due to the 
fact that they are calculated from the base functions C that 
are not error-free. Calculating the error propagation through 
a matrix inversion is extremely cumbersome. Furthermore 
the convolution in Eq. (3) implies that the error estimate 
should involve the integral of the unknown TF even if the 
errors in the individual measurements of the continuum flux 
are uncorrelated. Finally, adding error terms in the defini- 
tion (10) for the cross-correlation matrix of the base func- 
tions C means that the matrix E will not be a diagonal ma- 
trix. This is not a problem for the algorithm but the error 
estimate is likely to be sensitive to the interpolation scheme. 

There is a simpler way to measure the uncertainty due 
to the measurement errors in the continuum flux. Once the 
coefficients c have been obtained they are regarded as known 
with an arbitrarily high precision. The coefficients are com- 
bined with the line-flux measurements in the usual way and 



the error estimate is obtained from (13). The new element 
in the analysis is to combine the coefficients with the errors 
in the base functions C. This will produce an error estimate 
for each point of the integration kernel IC, with the propa- 
gated errors obtained from the interpolation scheme and the 
coefficients. Calculating this error propagation is relatively 
simple since it results entirely from linear combinations with 
known coefficients of the least-squares fits of the interpolat- 
ing polynomials which have associated error estimates (cf. 
Eq. 4). 

The departure of IC from the target form T is already 
used for an estimate of the error due to an imperfect match 
of the two (cf. Pijpers & Thompson, 1994). Instead of a 
strict upper limit to this error that would be obtained if 
the continuum flux was error free, the upper limit is now 
'fuzzy'. This arises because there is an error 5IC{t) associated 
with the integration kernel K,. Because the kernel is only 
used in integrated form the effect of a variance is reduced. 
If by analogy to the definition used by Pijpers & Thompson 
(1994) a 'systematic error' measure x is defined 

X = J dr [IC-T]\ (14) 



the upper limit to the systematic error in is given by : 

Kpp = y j (15) 

where 'I'min and ^max refer to the minimum and maximum 
value of the TF over the range [0, Tmax] ■ Only if this error 
i'upp is of the same order of magnitude as the a (/C(t)) or 
smaller than that, will the errors in the integration kernel 
significantly affect the error estimate -Eupp due to imperfect 
matching of the calculated integration kernel IC to the target 
function T. 

In this way the errors in the data are fully taken into 
account. The errors in the line fluxes are propagated prop- 
erly into an error estimate for the TF determination. The 
errors in the continuum fluxes arc propagated properly into 
an error estimate of the integration kernel IC. The depar- 
ture of K, from T can be used to estimate the effect of this 
on the estimate of the TF at the appropriate tq. It should 
be noted that this measure x is not available in the classi- 
cal MOLA methods such as the original Backus & Gilbert 
method (1967, 1968, 1970). This is due to the fact that the 
concept of a target function is unique to the SOLA method. 

At this point one could argue that the minimization in 
Eq. (7) is not optimal if ordy the errors in the line flux are 
included in the matrix E. From Eq. (3) it can be seen that if 
there is an uncorrelated measurement error associated with 
the continuum measurement as well as the flux measure- 
ment, then this will introduce an additional variance of the 
order of ^^a^ {C{ti)). The integral $ of Vl/ can be determined 
independently (cf. sec. 4.5). The combination of the errors 
in the continuum and emission-line light curves can be es- 
timated very roughly by adding the two. For the matrix E 
used in the optimization such a rough estimate is sufficient. 
The form that is used here is given by : 

= [a\L{U)) + '^'a\C{U))\/E^ 

Eij = i^j 

Because of the interpolation of the continuum light curve the 
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errors 5C are not really uncorrelated. It is therefore not cor- 
rect to use these estimates for the final error determination. 
It is therefore never used for that but instead the analysis is 
used that is described in the first part of this section. How- 
ever, for the minimization in (7) it is a sufficiently accurate 
approximation to ensure that the errors on the constructed 
kernel IC remain sufficiently small to produce meaningful re- 
sults. 

In Eq. (16) the factor En is introduced to produce a 
dimcnsionless quantity. In helioseismology it is customary 
to make the trace of the matrix E equal to the number of 
base functions. The same is done here : 

N 
i=l 

3.6 Simultaneity of line and continuum 

measurements 

In the previous sections it is assumed that the measurements 
of the line flux and of the continuum flux are simultaneous. 
This is not necessary for the method to work. As was men- 
tioned before at the start of the time series the line-flux 
measurements arc of little use so they could be omitted. 

If there are continuum-flux mccisurements not accompa- 
nied by a measurement of the line flux these can simply be 
incorporated into the base functions C. This does not add a 
new base function but it does improve locally the accuracy 
of the interpolation. 

If there are line-flux measurements not accompanied 
by a measurement of the continuum flux these can also be 
added to the data. In this case a new base function is added 
with a value at r = of this base function that is inter- 
polated between subsequent and successive measurements 
of the continuum flux. It is important to note that adding 
such line-flux measurements can serve to improve the S/N 
ratio in the TF somewhat, but that they cannot improve the 
resolution with which the TF can be determined since the 
appropriate time resolution in the continuum flux is absent. 

4 SIMULATIONS 
4.1 The set-up 

To test the SOLA method on typical AGN variability data 
we chose a representative continuum light curve and a repre- 
sentative TF to simulate real data. The representative con- 
tinuum light curve was taten to be the NGC 5548 light curve 
as published by Peterson et al. (1994). This four-year light 
curve was interpolated linearly onto a regular grid and then 
smoothed with a moving-average filter of total width 40 d. 
The weights for the individual points were taken to be equal. 
Note that this smoothing has implications for the time reso- 
lution that can then be attained in the determination of the 
TF. 

Previous MEM reconstructions of the TF yielded TFs 
with little or no response at zero lag and peaking at 10-20 d 
lag. We thus decided to use a TF with similar characteristics. 
We calculated a TF using the RAN model from Blandford 
& McKee (1982) where we chose an inner radius ro of the 
BLR of 15 light days. Thus our artificial TF peaks at 30 d 
lag and has little response at zero lag. The maximum lag 



Figure 2. The TF that was used for the simulations 

Tmax with a non-zero TF is 60 d. This TF was normalized 
to unity in order to compare the integral over the TF of 
the reconstructions with unity. The resulting TF is shown 
in Fig. 2. The sharp peak may be unrealistic but does pro- 
vide us with an estimate of how well the SOLA method can 
reconstruct sharp features. The peak was chosen to lie at 
30 d instead of 20 d because the smoothing of the contin- 
uum time series degrades the final attainable resolution and 
we decided to 'stretch' the simulated TF by 50% in order to 
retain approximately the same resolution relative to Tmax as 
in the real data where Tmax ~ 40 d. 

We then convolved the artificial continuum light curve 
with this TF and obtained an artificial emission-line light 
curve. The two light curves were then chopped such that 
only the second year of data (cf. Peterson et al. 1994) re- 
mained. This reduced data set was then resampled irreg- 
ularly onto the same grid as the original observations were 
done. In this way two ideal (i.e., noise- free) light curves were 
produced. 

To examine the effects of data errors, random noise from 
a Gaussian distribution with a standard deviation of 2%, 
5% and 10% of the signal was added. The routine RANI 
from the Numerical Recipes (Press et al. 1992) was used as 
the random-number generator and for every continuum light 
curve the same seed was used, such that the noise pattern 
was identical for each continuum curve except for the mag- 
nitude of the noise. The same was done for the emission-line 
light curves albeit with a different seed to avoid correlated 
errors between the continuum and emission-line time series. 
The eight light curves (four continuum and four emission- 
line with 0%, 2%, 5% and 10% noise each) are presented in 
Fig. 3. 

These light curves were then used to test the SOLA 
method. Ten simulation runs were devised according to the 
scheme in Table 1. In simulations 1-4 the noise was increased 
from to 10% on both light curves. These simulations reflect 
real data. In simulations 5-7 only the noise in the emission- 
line light curve was increased from 2 to 10% in order to 
test the effect of the errors on the emission-line data only, 
and in simulations 8-10 only the noise in the continuum 
was increased from 2 to 10% whereas the emission-line light 
curve was taken to be error-free. These last simulations thus 
measure the effect of continuum-error propagation. 

4.2 The results 

We have run the simulations twice, once with a high resolu- 
tion of A = 5 d and once with a resolution of A = 10 d. The 
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The eight light curves used for the simulations. The left column shows the continuum light curve and the right column the emission- 
line light curve. The upper row shows the noise-free curves, the second row 2% noise has been added to the data, the third row 5% noise, 
and the lower row 10% noise 



moan sampling interval of the time scries is only 3.4 d and 
one might argue for being able to reach this resolution but 
as explained in Sect. 3 the introduction of noise prohibits 
this. For both simulation runs, the trade-off parameter fj, 
between the importance of resolution on the one side and 
error suppression on the other side was set to 0.001. The 
results of the SOLA inversion on the simulated data are 
presented in Figs. 4 and 5. Figure 4 shows the 5-day reso- 



lution solution whereas Fig. 5 shows the 10-day resolution 
solution. As an example, Fig. 6 displays the kernels of SIM 2 
(i.e., the simulation with 2% noise on both the continuum 
and emission-line light curves) with A = 10 d. 

We used Nb = 2 and the expected maximum resolution 
for the error-free data due to the irregularity of the time 
series and the interpolation scheme used to compensate for 
that is Amin ~ 10 d (cf. Eq. 12). We thus expect the simu- 
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Figure 4. The recovered TFs with a five-day resolution for the simulated data. The target TF is shown as a dashed line in each plot 



lation with A = 5 d to be undersampled. 

The first thing wc note from Figs. 4 and 5 is that the 
error bars on the A = 5 simulation are systematically larger 
than on the A = 10 simulation. This is an obvious result 
due to the fact that less light-curve data are used per TF 
resolution element because the kernels K. for the A = 5 
simulation are narrower than for the A = 10 simulation. 

The second thing we note in both (A = 5 and A = 10) 
ideal simulations (SIM 1) is that the TF is recovered to 



within the error bars but shows a slightly displaced peak 

towards shorter r. This is due to the asymmetry of the used 
TF and the skewness of the kernels at r 20 d. The skew- 
ncss of the kernel can be seen from Fig. 6 where the kernels 
of SIM 2 with A = 10 d are displayed. The kernels that 
are significantly affected, at r = 15 d and 20 d, are no- 
tably asymmetric: they have more power towards larger lag 
than shorter lag. Together with the asymmetric simulated 
TF used it is obvious that this would raise the reconstructed 
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Table 1. The simulations 



continuum noise 


Emission-line noise 






0% 


2% 


5% 


10% 


0% 


SIM 1 


SIM 5 


SIM 6 


SIM 7 


2% 


SIM 8 


SIM 2 






5% 


SIM 9 




SIM 3 




10% 


SIM 10 






SIM 4 



TF above the input TF at these lags. 

The T = boundary of the reconstructed TF is slightly 
raised compared to the input TF, especially in the A = 10 
case, because there is no data on the negative side of the 
boundary and the finite width of the kernels let more- 
positive signal from larger t leak into the t = domain. 
If the peak of the TF would lie close to a resolution element 
away from r — 0, this boundary effect will pull the observed 
peak of the TF towards t = 0, resulting in an underesti- 
mation of the size of the BLR. The effect of the boundary 
is made clear by the kernels in Fig. 6. Here it is observed 
that with A = 10 d the kernel at t = 5 d is still effectively 
a kernel at r = d and also the kernel at r = 10 d is 
significantly affected by the boundary. This boundary effect 
must be taken into account when examining short lags in 
reconstructed TFs. 

Having established these artificial efi'ects on the ideal 
data set (SIM 1) we can increase the noise on both con- 
tinuum and emission-line light curves to 2%, 5% and 10% 
via simulations 2-4. The reconstructed TFs show increas- 
ing error bars but also substructure that must be entirely 
due to noise. In the undersampled A = 5 simulations we 
would hardly recognize the TF in the 5%-noise case (SIM 3) 
even though the reconstruction is consistent with the in- 
put TF to within the errors. SIMS is still acceptable in the 
A = 10 case though also here we undersample because the 
noise degrades the maximum attainable resolution of 10 d 
somewhat. The real data on NGC5548 have error bars on 
the 3% level, close to SIM 2 and significaiitly better than 
SIM 3. For the real data, a resolution of A ~ 10 d would be 
expected with A^s = 2. Even though the sampling rate of the 
NGC 5548 light curves is significantly better than 10 d it is 
not legitimate to aim to resolve structure in the TF on scales 
less than ~10 d! Actual experimenting with A^s and A on 
the real data suggested A^s = 2 and A = 8 is the preferred 
parameter setting for the NGC 5548 data. 

Simulations 5-7 show the effect of increasing the noise 
on the emission-line light curve while using the error-free 
continuum data. It is clear that this will increase the error 
bars and substructure considerably and it merely shows that 
the resolving power of the inversion degrades, as expected, 
and we aim for a too high resolution when we set A = 10 d. 
The error bars in SIM 7 are larger than in SIM 4 because of 
the way the data errors are used in the minimization. In the 
minimization the variance of the continuum measurements 
and line measurements are added. The error bars in the TF 
are calculated on the basis of the line-flux measurements 
only whereas the errors in the continuum flux are used to 
obtain an error estimate for the integration kernel. So with 
the same weighting factor n the trade-off in the minimization 
shifts towards smaller error variance in the TF from SIM 7 



to SIM 4. The error variance in the constructed integration 
kernel K is larger in SIM 4 than in SIM 7 (cf. Fig. 7) which 
compensates for the decrease of a (TF) . 

The simulations that used error-free omission-line light 
curves but increasingly noisy continuum light curves (simu- 
lations 8-10) show the eflFect of the error propagation 'under 
the integral'. These errors clearly introduce spurious signals 
and structure in the TF on seemingly random time scales 
larger than the resolution. This is a problem that one must 
test for in order to interpret a noisy inversion correctly. The 
effect, however, is moderate up to the 2-3% noise level. 

4.3 The resolution of the TF 

Figure 7 presents all kernels peaking at r = 30 d for both the 
A = 5 (left column) and A = 10 simulations (right column). 
These are the kernels responsible for picking up the TF sig- 
nal at T = 30 d. Obviously, these kernels have finite width 
and the full width at half maximum (FWHM) is a measure 
of the resolution of the reconstructed TF. Note that the 
FWHM is at best 2Av1n2. Equation (12) shows that with 
the present interpolation scheme, the maximum resolution 
attainable is Amin ~ 10 d for the present data set, even 
in the absence of noise! The FWHM of the kernel for the 
A = 5 SIM 1 simulation is ~ 30 d, implying Amin ~ 18 d. 
This is due to the fact that we cannot resolve high-frequency 
structure that is not present; recall that these simulations 
arc smoothed versions of a real light curve with a smooth- 
ing window of 40 days. For ideal data, the interpolation can 
be done better than done at present because A's can be de- 
creased. 

The kernels clearly show that increasing the noise on the 
continuum-fiux data increases the uncertainty in the deter- 
mination of the kernel, whereas the emission-line noise has 
no effect on the determination of the kernels. When increas- 
ing the noise, the kernels for the A = 5 simulations become 
narrower than the theoretical limit of Amin- The reason for 
this is that the noise is then used by the method just as real 
rapid variations would be used : to produce a small width 
of the kernel. The noise then allows a closer resemblance of 
the kernel to the target function. Thus, in the presence of 
noise, aiming for a resolution better than Amin means risk- 
ing spurious structure to show up in the final TF, because a 
narrow kernel also means fewer (noisy!) data points are used 
for the inversion. A spurious datum can thus have a large 
impact on the inversion when aiming for too high a reso- 
lution. This is effectively equivalent to undersampling the 
data. Comparing the A = 5 simulations with the A = 10 
simulations clearly shows this undersampling problem is not 
severe if we confine ourselves to A --^ Amin ~ 10 d. 

Note that aiming for a A less than Amin does not in- 
crease the error estimate of the TF. The error bars on the TF 
are calculated from the emission-line profiles only, whereas 
-Eupp is calculated from the mismatch of the kernel and the 
target function, in the case of noisy data, the noise allows 
the kernel to closely resemble the target function even in the 
absence of intrinsic high-frequency variations (cf. Fig. 7), 
thereby keeping Eupp small. With noise-free data, it is not 
possible to construct a kernel that is narrower than the mini- 
mum variability time scale of the intrinsic variations. Aiming 
for A < Amin thus equals risking to fit noise instead of sig- 
nal; this leads to the spurious substructure in the TF, and 
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Figure 5. The recovered TFs with a ten-day resolution for the simulated data. The target TF is shown as a dashed line in each plot 



this is evidently not a wise thing to do. 

4.4 The cross-correlation function 

It is well-known that the position of the peak of the cross- 
correlation function (CCF) of the continuum and emission- 
line light curves is a measure of the size of the BLR (Gaskell 
& Sparke 1986; Gaskell & Peterson 1987; Edelson & Krolik 
1988; Robinson & Perez 1990). It is not difficult to work out 



(e.g., Bracewell 1986; Penston 1991) that the CCF is the 

convolution of the TF with the auto- correlation function 
(ACF) of the continuum time series. Because the ACF is a 
symmetric function, asymmetric properties of the TF, like 
e.g. the position of the peak, must be reflected in a similar 
way in the CCF. Of course, the reverse is also true. The 
error analysis described in the previous sections, can thus 
also be applied, at least qualitatively, to the CCF. 

Table 2 presents the peak values and the peak positions 
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Figure 6. The kernels for simulation 2 with a resolution A = 10 d 



of the CCFs of the time series used for SIM 1-SIM 10. It 

is clear that increasing the noise on either the continuum 
or emission-hue light curve decreases the correlation coeffi- 
cient of the peak of the CCF; noise increases the width of 
the ACF and the CCF peak becomes less well-defined. This 
is confirmed by Table 2. Also, noisy data decrease the corre- 
lation coefficient at small delays (provided the noise on the 
continuum and emission-line light curves are not positively 
correlated), thereby increasing the position of the peak of 



the CCF towards larger delays. This is also confirmed by 
Table 2. Thus, noisy data most likely overestimate the size 
of the BLR when only the peak of the CCF is considered. 
Note that theoretical modelling of BLRs show the peak of 
the CCF is more sensitive to the inner boundary of the BLR 
than to its centre, thereby underestimating the actual size of 
the BLR (Robinson & Perez 1990). However, this is different 
from the numerical effect presented here. 

Because negative lags are unphysical, with finite resolu- 
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Figure 7. The kernels for all simulations for r = 30 d. The left column displays the 5-day resolution case; the right column the 10-day 
resolution case. An artificial vertical offset has been given for displaying purposes 



tion not only the TF, but also the CCF will show a bound- 
ary eflFect at r = 0. The correlation coefficient at r = will 
thus be increased for TFs peaking away from zero delay, and 
this effect will be larger for worse resolutions. This pulls the 
peaks of the TF and CCF towards shorter lags. In the data 
set examined hero this might be a problem, especially with 
the A = 10 resolution, where the kernels are significantly 
asymmetric up to r ~ 20 d due to the boundary. Also, the 
asymmetric TF together with the finite resolution will shift 
the peak towards smaller lag. The peak of the TF is known 
to be at r = 30 d but both the TF and the CCF peak at 
T ~ 23 d for the noise-free data due to the low resolution of 
the experiment. 



Noisy and incomplete data thus reduce the resolving 
power of both the TF and the CCF such that the peak 
positions of the TF and CCF are displaced. 

4.5 The integral of the TF 

As Wanders (1994) showed, the delay integral of the TF is a 
physically relevant parameter $ = J dr\E'(r) that measures 
the total observed line reprocessing in the BLR. $ is not 
directly reproducible by the MEM; only by first solving for 
the TF and then integrating. 

From Eq. (6) it can be seen that an estimate $ of the 
integral $ can be obtained by letting the kernel IC be unity 
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Table 2. The cross-correlation results and the 
integral <J> of the TF for the simulations. The 
error in the peak position of the CCF is about 
2 d 















CCF(rpeak) 








SIM 1 


0.91 


23 d 


1.000 


0.1% 


SIM 2 


0.89 


23 


0.998 


0.3% 


SIM 3 


0.79 


25 


0.996 


0.7% 


SIM 4 


0.62 


37 


1.002 


1.4% 


SIM 5 


0.90 


23 


1.000 


0.3% 


SIM 6 


0.85 


29 


1.000 


0.7% 


SIM 7 


0.74 


34 


1.000 


1.4% 


SIM 8 


0.90 


23 


0.997 


0.1% 


SIM 9 


0.85 


30 


0.993 


0.1% 


SIM 10 


0.73 


30 


0.993 


0.1% 



over T € [0, Tmax]- This can be accomplished by aiming for 
an infinitely low resolution, i.e., by setting A equal to a very 
large number. In this way, SOLA optimizes correctly for the 
integral problem. The calculation can be done at any to with 
this fiat kernel. Similarly to the calculation of ^(to), a di- 
rect error estimate for $ is obtained. The results for the ten 
simulations are presented in Table 2. Except for the unreal- 
istic case of orror-frcc line flux measurements the deviation 
from 1 of the integral is smaller than the error estimate from 
SOLA based solely on the line-flux errors. This means that 
the errors in the continuum flux do not prevent an accurate 
measurement of the integral. It is immediately clear that $ 
can be recovered to a precision better than ~1% for realistic 
data with ^5% noise. For the real data on NGC 5548, which 
have characteristics close to SIM 2, we can therefore expect 
an accuracy of 3> of 

Constant contributions to the fluxes can confuse this 
result but will not be discussed here because they do not 
influence the conclusions of the present paper. See Wanders 
(1994) for a discussion of this. 

5 NGC 5548 

Peterson et al. (1994) published four years of monitoring 
of NGC 5548 and presented MEM-recovered TFs for each 
year of monitoring. We have taken the same data set and 
recovered the TF of each year using the SOLA method. We 
did not aim for a higher resolution than A = 8 d since 
this would introduce local structure in the TF that is due to 
measurement errors in the continuum and emission-line light 
curves. The trade-off parameter jj, was set to 0.001 for all four 
inversions, and the maximum r for which we calculated the 
TF was 50 d. Runs with larger Tmax showed no significant 
structure at larger lags. Since we have to 'throw away' the 
first data of the the emission- line light curve up to Tmax, 
decreasing Tmax means including more data. 

The results for the four individual years are shown in 
Fig. 8. These TFs resemble the TFs derived by the MEM 
(Peterson et al. 1994); little response at zero lag, little re- 
sponse at lags larger than about 40 d and peaking at approx. 
12-20 d. The decreasing trend in the position of the peak 
of the TF from ~ 20 d in 1989 to ~ 13 d in 1992 as was 
reported by Peterson et al. is also reproduced here (cf. Fig. 
8). The peak of the TF is always more than a resolution 



Table 3. The integral * of the TF 

for NGC 5548 



year 


'I'SOLA 


5*SOLA 


*av 


1989 


0.863 


0.004 


0.876 


1990 


0.797 


0.004 


0.825 


1991 


0.776 


0.007 


0.794 


1992 


0.720 


0.006 


0.738 



element away from the boundary and therefore the bound- 
ary will not have too much influence on the position of the 
peak. However, the simulations have shown that a system- 
atic underestimation of the position of the peak of ~5-7 d 
cannot be excluded. The upper limit Eupp to the system- 
atic error due to the uncertainty in the kernels is signifi- 
cantly smaller than the statistical uncertainty propagated 
from the line measurements, indicating we do not run into 
serious trouble due to the uncertainties in the continuum 
measurements. Note that the TF is not significantly away 
from zero at zero lag in any of the reconstructions, except 
for the year 1989. There is no indication that the TF would 
be significantly negative anywhere. This does not exclude 
the TF to become negative on smaller time scales than the 
present resolution of about 10 d, but the positiveness of the 
TFs suggests the MEM does not run into serious trouble 
while using its positivity constraint, at least not with the 
resolutions attained with the present data. Negative values 
in the TF might arise e.g. due to nonlinear line responses 
(Sparke 1993). 

Note that the SOLA method is not hampered or aided 
by a positivity constraint and negative responses would have 
been recovered as long as a convolution of a gaussian with 
such a negative signature in the TF is larger than the noise 
level. Experimenting with the resolution (A) showed that 
when a too high resolution was aimed for, artificial structure 
appeared in the TFs, similar to the noisy simulations in 
Sect. 4. These artificial structures can become significantly 
negative, but are totally due to noise. Care must be taken 
not to overinterprct such structure. 

In all, the present SOLA results are consistent with 
previous MEM results for the recovery of the shape of the 
NGC 5548 TF. This gives us more confidence in the viabil- 
ity of these results and it presents the SOLA method as an 
alternative to the MEM. 

Table 3 presents the integral of the four TFs of 
NGC 5548 as derived by the SOLA method ($sola) and 
as published by Wanders (1994) who derived <I> via the av- 
eraging method $av ~ {L{t))/{C{t)). The errors on <I>av are 
approximately 2.5% and as discussed in Sect. 4.5, the errors 
on $soLA are on the order of 1%. It can thus be seen that 
the decreasing trend in $av is confirmed by #sola. 

A decrease in $ is equivalent to a decrease in observed 
reprocessing of the emission-line within the BLR (Wanders 
1994). We thus confirm the findings that both the position 
of the peak of the TF and the total integral of the TF are 
time-variable on a time scale of a few years, corresponding 
to the dynamical time scale of the BLR clouds. This strongly 
suggests evolution of the BLR takes place in NGC 5548 and 
the system is not stationary on these time scales. 
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Figure 8. The four transfer functions for the H/3 emission-hne region of NGC 5548. The TFs are for the monitoring years 1989, 1990, 
1991 and 1992. Also shown are the upper limits Bupp to the systematic error made due to the uncertainty in the kernels (cf. Eq. 15) 



6 RESUME 

This paper is primarily concerned with the development 
and testing of a new method for the inversion of time-series 
of AGN contiimum and line-emission fluxes, to obtain the 
transfer function of the broad-Une region. The results of the 
application of this method to real data of NGC 5548 are also 
shown. The various aspects of the method, of the tests and 
of the apphcations are discussed extensively already in sec- 
tions 3, 4, and 5 respectively so we will confine ourselves to 



a few summarizing statements. 

The SOLA method constructs an integration kernel 
from a linear combination of subsections of the continuum 
light-curve. The TF is reconstructed from a linear combina- 
tion of the emission-line measurements corresponding to the 
continuum subsections, using the same linear coefficients. 
The value of the TF thus obtained is a weighted average 
of the TF where the weighting function is the constructed 
integration kernel. 

The method is not subject to any regularity constraint 
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or positivity constraint for the TF. 

The method allows a very clear determination of the 
time resolution with which the transfer function can be de- 
termined. Most of the procedure depends only on the sam- 
pling strategy and the standard deviations in the flux mea- 
surements. It is therefore possible to determine the resolu- 
tion before reconstructing the TF. This provides a-priori a 
criterion for avoiding or discarding spurious sub-structure in 
the TF. 

The method uses the variance in the data to produce the 
smallest possible standard deviation in the TF at the time 
resolution used. Because the method is linear the calculation 
of error propagation is particularly straightforward. 

Because only one large- matrix inversion has to bo done, 
the SOLA method is cost-effective and a typical inversion of 
the NGC 5548 data (typically ~80 base functions) takes ap- 
proximately one hour of CPU time on a SUN Sparc-2 work- 
station. This allows for experimenting in parameter space 
by varying the control parameters /i and A. In the origi- 
nal Backus- Gilbert formulation the large number of matrix 
inversions prevents parameter-space exploration. 

The tests of the method show that for a realistic TF and 
time sampling the TF can be recovered to within the formal 
error estimates. Because of the finite resolution of the time 
series some systematic shifts axe found in the position of the 
peak of the TF. We argue that this is essentially an effect 
of resolution. It is therefore intrinsic to the problem and the 
sampling. This means that other methods to determine this 
peak position such as the CCF and the MEM must also 
suffer from these systematic effects. 

The application of the method to the data obtained for 
NGC 5548 confirms the time evolution found previously of 
the TF over the four years 1989-1992: both the position of 
the peak of the TF and its total integral $ are time- variable. 
This evolution of the broad-line region takes place on the 
dynamical time scale of the clouds within the BLR. 

The results of the method, the tests, and the applica- 
tion to real data allow us to make precise recommendations 
for setting up observing campaigns, based on the required 
resolution of the TF and the expected standard deviation 
in the flux measurements. The strategy of current observing 
campaigns cannot be oxpoctod to produce rrmch higher reso- 
lution or better results than shown here and we strongly urge 
future monitoring campaigns to use a higher sampling rate 
and an even better photometric accuracy. Unfortunately, 
both will increase the pressure on telescope time. 

The method can easily be adapted to a 2-D inversion for 
the TF as a function of time-lag and velocity. Considering 
that the shape of an emission-line profile can be determined 
more accurately than its flux this may yield a great deal of 
new information. For instance, the SOLA method opens the 
way to the determination of the Q-function as introduced 
by Perry, van Groningen & Wanders (1994). They split the 
full TF into a steady-state and a deviation term: 

•^{v,T)^^^{r) + Q{v,T) (18) 

where $(v) = / dT*(v,r) is the steady-state profile, i.e., 
the emission-line profile in response to a constant continuum 
flux. Q{v,t) is then a function extremely sensitive to varia- 
tions in the emission-line profile and from this a great deal 
about the structure of the BLR can be learned. Perry et al. 



present an inversion problem to directly measure Q{v, t) via 
normalized emission-line profiles. Because Q{v, r) is partly 
negative, the MEM cannot recover this function because of 
its built-in positivity constraint. This is no problem for the 
present SOLA method. 

Also the steady-state profile $(«) is an important quan- 
tity and can be measured from variable sources by the SOLA 
method similarly to the way $ is measured as described in 
Sect. 4.5. Evolution of the BLR wiU reflect itself in <!>(«) 
also and $(«) will thus be a powerful way to study evolution 
while simultaneously being comparatively easy to calculate 
by the SOLA method. 

The SOLA method is thus not only an alternative to 
existing inversion techniques for the derivation of the TF, 
but can be expected to be a powerful tool in the study of 
emission-line profiles as well, which opens the way to a bet- 
ter understanding of the broad-line region in active galactic 
nuclei. 
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